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Abstract 

We discuss the interpretation of the annual modulation signal seen in the DAMA experiment in 
terms of spin-independent elastic WIMP scattering. Taking into account channeling in the crystal 
as well as the spectral signature of the modulation signal we find that the low-mass WIMP region 
consistent with DAMA data is confined to WIMP masses close to ~ 12 GeV, in disagreement 
with the constraints from CDMS and XENON. We conclude that even if channeling is taken 
into account this interpretation of the DAMA modulation signal is disfavoured. There are no 
overlap regions in the parameter space at 90% CL and a consistency test gives the probability 
of 1.2 X 10^^. We study the robustness of this result with respect to variations of the WIMP 
velocity distribution in our galaxy, by changing various parameters of the distribution function, 
and by using the results of a realistic A^-body dark matter simulation. We find that only by making 
rather extreme assumptions regarding halo properties can we obtain agreement between DAMA 
and CDMS/XENON. 
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I. INTRODUCTION 



The DAMA collaboration has collected an impressive amount of data in their search for 
the scattering of weakly interacting dark matter particles (WIMPs) off Sodium Iodine. The 
combined data from DAMA/Nal (7 annual cycles) and DAMA/LIBRA (4 annual cycles) 
amounts to a total exposure of 0.82 ton yr in a field where exposure is measured in units 
of kg days. DAMA/LIBRA has now provided further evidence for an annual modulation of 
the event rate in the energy range between 2 and 6 keVee, the claimed statistical confidence 
of the positive signal being 8.2a [|T|. The phase of the observed modulation (with maximum 
on day 144 ± 8) is in striking agreement with the expectation for a modulation in a WIMP 
scattering signal due to the rotation of the Earth around the Sun, (expected maximum 
day 152, June 2nd), see e.g., |^ for a review. An interpretation of this effect in terms 
of spin-independent interactions of conventional WIMPs with masses > 50 GeV is 
in direct conflict with the constraints from several experiments looking for direct WIMP 
detection, most notably with the data from CDMS and XENONIO which exclude 
the WIMP cross section consistent with the DAMA modulation for ~ 50 GeV by many 
orders of magnitude. In light of this, several alternative explanations of the DAMA annual 
modulation have been proposed, for example spin-dependent interactions |P, light WIMPs 
with < 10 GeV masses 0, ||, 0, keV scale axion-like dark matter |jlO[ (see however, [jlT], 0), 



dark matter interacting only with electrons [r^|, inelastic WIMP scattering [0, and 
mirror dark matter |jl6| . 

In this work we reconsider the possibility of spin-independent elastic scattering of light 
WIMPs with < 10 GeV masses 0, |, §, see [|T7|, 0, 0, |^ for recent studies. The original 
idea is that light dark matter scattering on the relatively light Sodium nuclei in DAMA 
could deposit enough energy in the detector to give a signal, whereas the scattering of light 
halo particles off heavier nuclei, such as for example Ge in CDMS or Xe in XENON would 
lead to energy depositions below the threshold of those detectors. Recently the importance 
of the so-called channeling effect pT| in the crystal structure of the experiment has also been 



emphasized ITR Specific models for WIMPs with 
example in O, E^, E^, 



10 GeV have been studied for 



J^, 3 . Here we do not discuss theoretical implications but focus on the 
phenomenology of direct detection experiments in a model-independent way by assuming 
that such light WIMPs can provide the correct relic abundance while any direct collider 
constraints can be evaded. 

In this region of WIMP masses several experiments p5| , P^ , ^ exclude WIMP-nucleon 
scattering cross sections in the range ap > 10~^° cm^. As we will see in the next pages, 
once we have included channeling as well as the spectral shape of the DAMA modulation 
signal, the allowed region of our interest is obtained at much small cross sections, around 



10 cm^ and m 



^ 10 GeV. In this region the most relevant constraints come 
the 2008 Germanium data from CDMS HI, and the 2005 CDMS data on 



from XENON _ _ 
Silicon [E^ . Indeed, as we will discuss, the spectral shape of the DAMA annual modulation 



restricts and ap to a region excluded by these experiments. 

In our study we elaborate on this result and discuss how robust it is with respect to 
different assumptions about the dark matter halo of our galaxy. The impact of non-standard 
halo properties on dark matter direct detection experiments has been discussed by many 
authors, see for example [BO, 31, 32, RSl. At a qualitative level, one would expect that smaller 



2 



velocity dispersions or truncated velocity distributions would seem to favour the dark matter 
interpretation of the DAMA signal, as they could lead to more events above the low energy 
threshold of DAMA but below that of other experiments. Furthermore, anisotropies in the 
velocity dispersion could amplify annual modulation signals. 

The outline of our work is as follows. In Sec. ^ we briefly summarise the phenomenology 
of elastic WIMP scattering in direct detection experiments and give some technical details 
on our analysis of DAMA, CDMS and XENON data. The results for a standard dark 
matter halo are presented in Sec. |T|. In Sec. we consider deviations from the standard 
assumptions made about the WIMP velocity distribution: we use results from the Via 
Lactea A^-body dark matter simulation [Q, we vary several parameters of the Maxwellian 
distribution and consider asymmetric velocity profiles. Sec. ^ contains our conclusions. In 
Appendix ^ we comment on the DAMA fit using the annual modulation energy spectrum, 
and in Appendix ^ we briefly compare our results to the ones from other authors. 

II. THE WIMP SIGNAL IN DIRECT DETECTION EXPERIMENTS 

In this section we briefly summarise the phenomenology of WIMP scattering and describe 
our analysis of DAMA, CDMS and XENON data. 



A. The event spectrum from elastic WIMP scattering 

The differential event spectrum for WIMP scattering in counts per unit mass of a given 
nucleus per unit exposure time and per unit energy as a function of the recoil energy En is 
given by the expression (see e.g., p|) 

Here p is the local WIMP energy density for which we adopt the canonical value p = 
0.3 GeV/cm^, ap is the WIMP scattering cross section on a proton^, A is the mass number 
of the target nucleus, pp = m^mp/{m^ + mp) is the reduced WIMP-proton mass and we use 
the common Helm form factor F{q) = 3e~'^ * /^[sin(gr) — gr cos(gr)]/(gr)^, with s = 1 fm, 
r = a/-R^ — Ss^, R = 1.2A^/^ fm, q = \/2MEr, with M being the nucleus mass. The 
function r] contains the integral over the WIMP velocity distribution: 

//•oo 
dn^ dvvf^{w,t), (2) 



where fmin = a/ MEr/2p\j is the minimum velocity of a WIMP to produce a recoil energy 
Er, and V = |v|. The WIMP velocity distribution in the Earth rest frame /®(v, t) is obtained 
from the distribution in the galactic rest frame /gai(v) by 

/e(v,t) = /gai(v + V0 + Ve(t)). (3) 



^ Note that only the product of p x Cp is relevant for the scattering rate. Therefore, whenever we use the 
symbol ap the cross section is implicitly normalised to the value of p — 0.3 GeV/cm"^. 
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In the coordinate system in which x points towards the galactic center, y towards the 
direction of galactic rotation, and z towards the galactic north pole, we use for the velocity 
of the Sun = (0, 220, 0) + (10, 13, 7) km/s (including the local Keplerian velocity of 
220 km/s ||3^ as well as the Sun's peculiar velocity, see also and references therein). To 



describe the motion of the Earth around the Sun we use the parametrisation of |^: ^®{t) = 
%(eisinA - egcosA), with = 27rA.U./yr = 29.8 km/s, ei = (-0.0670,0.4927,-0.8676), 
62 = (-0.9931, -0.1170, 0.01032), and A(t) = 2n{t - 0.218). 

The "standard halo model" assumes for the DM distribution an isotropic isothermal 
sphere, which leads to a Maxwellian velocity distribution in the galactic frame, truncated 
at the escape velocity Vesc- 

= 1 V> .esc ' 

where we adopt as default values v = 220 km/s and fesc = 650 km/s. Here and throughout 
the paper we use the notation = 2(< v"^ > — < v >^) = 2cr^. In order to properly take 
into account the impact of the finite escape velocity as well as allowing for non-standard 
halos deviating from Eq. ^ we perform the integral in Eq. ^ numerically. In Sec. Ill we first 



consider the standard halo model, whereas in Sec. fV] we go beyond these default assumptions 
by varying the parameters of the velocity distribution as well as changing its shape. 



B. On quenching and channeling 



In the analysis of DAMA data the effects of quenching and channeling are important |21 



38| . For quenched events the recoiling nucleus loses its energy both electromagnetically as 
well as via nuclear force interactions, where the light yield in the scintillator comes mainly 
from the electromagnetic part. To take this effect into account the event energy is measured 
in equivalent electron energy (in keVee), defined by g x E^i for the total nuclear recoil energy 
in keV. For the elements in DAMA one has gNa = 0.3 and qi = 0.09. However, due to the 
crystalline structure of the target, for certain angles and energies of the particles no nuclear 
force interactions happen and the entire energy is lost electromagnetically. Hence, for these 
so-called channeled events one has g ~ 1, see For the fraction / of channeled events 



relevant for DAMA we use the parameterisation 

" 1 + 0.75E, ' ^^^^^^ " 1 + 0.65E, 

for En in keV. These expressions reproduce to good accuracy the curves shown in fig- 
ure 4 of Departing from Eq. |^, the predicted spectrum in DAMA (in units of 
counts/kg/day/keVee) is obtained by 

Rdama{E)= Y1 m Im {[^-fx{E/qx)]Rx{E/qx) + fx{E)Rx{E)} , (6) 

where the first term in the curled bracket corresponds to quenched events and the second 
to channeled (and therefore unquenched) events. 
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Channeling does not occur in liquid Nobel gases like in the XENON experiment. Since no 
information on channeling in Germanium and Silicon is available for us, we do not take into 
account channeling in CDMS. Note, however, that CDMS requires the coincidence of signals 
in phonons and ionisation and hence, since channeled events would not give a phonon signal 
they would not look like a WIMP signal defined by the coincidence. Therefore, the fraction 
of channeled events corresponds effectively to an efficiency factor reducing the effective 
exposure. Hence, if channeling was indeed relevant for CDMS the final exclusion limits 
would be somewhat weaker. 

In conclusion, channeling is an important effect for the interpretation of data from direct 
detection experiments and we stress the need of reliable information (probably requiring 
dedicated measurements) on this effect for any solid DM detector. 



C. Fitting DAM A/LIBRA data 

For the model-independent analysis of DAMA data the signal as a function of energy and 
time is parametrised as 

S{E,t) = So{E) + A{E) cos uj{t~ to), (7) 

with u = 27r/l yr, = 152 day. For our analysis we use the data on the modulation 
amplitude A{E) for the full 0.82 ton yr DAMA exposure^ given in figure 9 of in 36 bins 
from 2 to 20 keVee. As we will see the spectral shape of the signal is quite important for 
constraining the WIMP parameters. The prediction for the modulation amplitude in an 
energy bin i from E~ to E^ is obtained from Eq. ^ by 

Ar'' = J dE^[Rr,AMA{E, t = 152) - Rbama{E, t = 335)] ' dE'G{E, E') (8) 



where G{E, E') is a Gaussian energy resolution function with width []39 

^DAMA/^ = 0.45/ [keVee] + 0.0091 . (9) 
Then we construct a function 

xLmaK, = E (^ ^-"^""^^^'^-" ) ' , (10) 

using the experimental data points Af^^ and their errors o"j from figure 9 of [Q.^ We find 
the best fit point for the WIMP mass and the scattering cross section by minimising Eq. |l^ 
with respect to and Up. Allowed regions in the (m-^, Up) plane at a given CL are obtained 



^ Here and in the following we use the acronym "DAMA" to denote the combined DAMA/Nal + 

DAMA/LIBRA data, except where explicitly noted otherwise. 
^ Fig. 10 of [Q shows that the Af^'^ are consistent with being Gaussian distributed, justifying the adopted 

in Eq. ^ 
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by looking for the contours x^i^x^ ^p) ~ Xmin + ^X^(CL), where Ax^(CL) is evaluated for 
2 degrees of freedom (dof), e.g., Ax^(90%) = 4.6 or Ax^(99.73%) = 11.8. 

In general the constant part of the spectrum, So{E), will consist of a time-averaged dark 
matter contribution (R) plus an un-identified background B: 

SoiE) = {R{E)) + B{E) . (11) 

In a given model such as for example WIMP scattering, the annual modulation amplitude 
A{E) and the averaged signal {R{E)) are not independent. Hence, for a given fit to the data 
on A{E), the expected constant signal {R{E)) can be predicted by using Eq. |]. In order 
to take this additional information into account we use the data from figure 1 of which 
shows the constant signal 5*0 in 32 energy bins from 2 to 10 keVee for the DAMA/LIBRA 
detectors. For each pair of {m^, ap) we calculate the expected signal from WIMP scattering 
(R) in each of these energy bins. Whenever (R) exceeds the observed rate in one of the bins 
that particular values of {m^, ap) are not consistent with the data and have to be excluded. 
Note that for event rates of order 1 count/kg/day/keVee and the DAMA/LIBRA exposure 
of 0.53 t yr statistical errors are negligible for this purpose. 



D. Analysis of CDMS and XENON 

In our analysis we include the constraints from CDMS 2005 data using Silicon (CDMS- 
Si) 1^ , which, despite the relatively low exposure of 12 kg day, provides good sensitivity to 
the low-mass WIMP region because of the light mass of the target nucleus (Msi ~ 26 GeV) 
and the low analysis threshold of 7 keV. Furthermore, we include CDMS 2008 data on 
Germanium (CDMS-Ge) with a threshold of 10 keV. We use the energy dependent 
efficiency from Fig. 2 of which reduces the total exposure of 398.7 kg day to an effective 
exposure of about 121.3 kg day. For both, CDMS-Si and CDMS-Ge, no event has been 
observed. We calculate the expected number of events N^'^'^'^ as a function of and ap 
by integrating Eq. ^ over the relevant energy range and scaling with the exposure. A is 



constructed using the common expression for Poisson distributed data [Q, which for zero 
observed events simply becomes 

X^DMS = 2iVP-<^ . (12) 

Exclusion contours are defined by the standard Ax^ cuts for 2 dof with respect to the 
minimum, which of course occurs for N^^'^'^ = 0. Conceptually this prescription differs from 
the usual way to set a limit on ap for fixed by requiring A^p'''"^ < 2.3 for a 90% CL limit. 
However, by accident, since Ax^(90%) = 4.6 for 2 dof, in practice our definition in Eq. ^ 
leads to the same exclusion contour as the more conventional method of setting a limit on 
ap. 

For the analysis of data from the XENONIO experiment (XENON for brevity) we 
proceed in the following way. Using the 7 bins in nuclear recoil energy from 4.5 to 26.9 keV 
of table 1 of [§] the predicted number of events in bin i, Nf'^'^'^{m^, ap) can be calculated by 
integrating Eq. || and scaling with the exposure 316 kg day as well as the bin dependent 
efficiencies ec and Anr given in table 1 of [ffl . After the publication of Ref . M the so called 



parameter C^s relevant for the nuclear recoil energy scale in XENON was remeasured pO 



Whereas in H] a constant value Ces = 0.19 was used. Fig. 3 of |40| shows an energy dependent 
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deviation of £cfr from that value. We use the information from this figure to correct nuclear 
recoil energies in XENON. This leads to a somewhat higher energy threshold of about 5.5 keV 
(instead of 4.5), which shifts the bound on DM parameters to slightly higher values of m^. 

XENON observes 10 candidate events whose recoil energies can be inferred from figure 3 
of 1^. They are distributed over the 7 bins as (A) = (1,0,0,0,3,2,4), with an expected 
background (5^) = (0.2,0.3,0.2,0.8,1.4,1.4,2.7). We use the for Poisson distributed 



data 106 



XXENON 



N 



pred 



5,- A + Afog 



(13) 



where the second term in the square bracket is zero if = 0. Again we define the exclusion 
curve in the {m^, ap) plane by Ax^(CL) contours for 2 dof with respect to the minimum. 

In both cases, CDMS and XENON we include an energy resolution of 20%/ a/ Er [ke V] , 
and an uncertainty in the energy scale of 10%, which is added to the definitions Eqs. ^ 
and |T3| with the help of nuisance parameters. Note that CDMS and XENON report their 
results directly in terms of the recoil energy already corrected by the quenching factor. In 
contrast to DAMA, here this is possible because only a single element is used as target. 



III. STANDARD HALO RESULTS 

Figure ^ summarises our results assuming standard halo properties, showing the allowed 
region from DAMA together with the constraints from CDMS-Si, CDMS-Ge and XENON. 
First we discuss the fit to DAMA data alone (without constraints from CDMS and XENON). 
We find two islands in the {m^, ap) plane where DAMA can be accommodated. The best fit 
point is obtained at 

= 12 GeV , = 1.3 x 10"^^ cm^ , XoAMA.min = 36.8/34 dof , (14) 

with an excellent goodness of fit of 34%. There is also a local minimum at = 51 GeV 
with xfocai — 47.9. This solution is disfavoured with respect to the best fit point at about 
3cr for 2 dof (Ax^ = 11.1). The allowed regions around ~ 50 GeV shown in figure |l| 
are defined with respect to the local minimum. The low and high WIMP-mass solutions 
correspond to channeled and quenched scatterings on Iodine, respectively. In contrast to 
the situation when all events are assumed to be quenched 0, it turns out that scattering 
on Sodium is not relevant once channeling of Iodine events takes place. The reason is that 
quenched events on Sodium require a similar WIMP mass as channeled events on Iodine 
(i.e., ~ 10 GeV) but a much larger cross section ap (due to the dependence of the 
total cross section on the nucleus), and therefore, are highly suppressed once channeled 
scattering on Iodine takes place. In principle there would be also a solution from channeled 
events on Na, around ~ 5 GeV. However, it turns out that in this case the un-channeled 
events on Na still contribute to the signal, and indeed prevent fitting the data with the 
channeled Na events. Note that the solution around 50 GeV is excluded by some 

orders of magnitude by XENON and CDMS-Ge, and therefore we focus in the following on 
the low-mass region ~ 10 GeV. 

The gray contours in figure |^ correspond to an alternative method of fitting DAMA. 
Instead of using the detailed spectral information of the annual modulation, we fit the time 
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FIG. 1: Allowed regions at 90% and 99.73% CL for WIMP mass and scattering cross section on nucleon 
for DAMA, and exclusion contours for CDMS-Si, CDMS-Ge and XENON at 90% CL. We also display the 
limit from CoGeNT extracted from figure 2 of ||2^. The global best fit for DAMA is marked with a star, the 
allowed region around ~ 50 GeV is defined with respect to the local minimum, which is marked with a 
dot. For DAMA we show the regions obtained from using only the modulation amplitude for 2-6 keVee (gray 
curves) and from using the spectral shape of the modulation signal (shaded regions) . For parameters above 
the dashed curve the predicted number of events in DAMA/LIBRA is larger than the observed number of 
events. 



dependence of their signal integrated over energy. In figure 6 of |T]] data on the residual 
rate {S(t) — Sq) (c.f., Eq. 0) is given in 7 time bins of one single annual cycle. For the gray 
contours in figure |1] we use these data for the energy intervals 2 to 6 keVee and 6 to 14 keVee, 
where in the latter interval data are consistent with no annual variation. These results are 
very similar to the ones of (if channeling is neglected, not shown in the figure) and |19 



(including channeling), where only two data points for the modulation amplitude below and 
above 6 keVee have been used. 

We observe from figure || that the two methods of analysing DAMA data are consistent 
with each other (as it should be), but also that using the spectral information gives sig- 
nificantly stronger constraints on the allowed region. This is illustrated in figure ^ (left), 
showing the 36 data points on the modulation amplitude Ai used in our default analysis. 
While the prediction from the best fit point of Eq. nicely follows the data (solid curve), 
moving to smaller WIMP masses leads to a modulation signal more peaked at the lowest 
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FIG. 2: Left: Energy distribution of the annual modulation amplitude from DAMA/Nal and 
DAM A/LIBRA data extracted from figure 9 of (points with error bars), together with the prediction 
for three examples of WIMP masses and scattering cross sections (curves). Right: Energy distribution of 
the time averaged rate observed in DAMA/LIBRA extracted from figure 1 of (points), together with 
the prediction for two examples of WIMP masses and scattering cross sections (thick curves) as well as 
the corresponding un-identified background (thin curves) . The data are corrected for the energy dependent 
efficiency. 

energies. Therefore, although it is still possible to obtain the integrated signal in the inter- 
val from 2 to 6 keVee, the spectral shape is clearly inconsistent with data, as illustrated for 
= 6 GeV by the dashed curve. 

Finally we mention the implication of the data on the time averaged rate observed in 
DAMA. Parameter values above the dashed curve in figure |I] are excluded because they would 
lead to a higher event rate than observed. This leads to additional constraints for the high- 
mass solution. In figure |^ (right) we show the observed rate together with the predictions 
for the two local minima. Note that for the DAMA/LIBRA exposure of 0.53 t yr statistical 
errors are not visible at the scale of the plot. Clearly, solutions predicting a relatively large 
rate require that the un-identified background drops rapidly in order to give space for the 
WIMP signal. In particular, the solution at = 51 GeV requires that the background 
drops to zero in the first energy bin. Although this cannot be excluded a priori, at least 
such a background shape seems somewhat unlikely. The issue is less severe for the best fit 
point at = 12 GeV, since the ratio of modulation amplitude to average rate increases for 
decreasing WIMP mass. However, any point close to the dashed line in figure |1| is affected 
by this problem. 



The value for the cross section dp = 2.8 x 10~^^ cm^ formally gives the best fit to the data shown in 
figure ^ (left) for = 6 GeV. However, the value required to obtain the integrated modulation amplitude 
for this is about a factor 2 larger, o-p = 6 x 10^"'^ cm^, as can be seen from the gray contours shown 
in figure ||. 
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From figure |T] we find that the parameters allowed by DAMA data at 90% CL are excluded 
by the 90% CL limits of CDMS-Si, CDMS-Ge, and XENON. If all data are combined by 
adding the individual functions, 



Xglobal 



— XdAMA + XCDMS-Gc + XCDMS-Si + XxENON ' (1^) 

we find the minimum at = 9.5 GeV and = 1.2 x 10^'^^ cm^ with Xgiobai,min — 59.3/ (45 — 
2) dof , which corresponds to a 5% goodness of fit. 

Let us note that the goodness of fit test based on x^jj^/dof often is not very sensitive to 
tensions in the fit, especially in case of a large number of data points. This can happen if 
there are many data points which actually are not sensitive to the relevant parameters, and 



hence, allow to "hide" the problem in the fit, see for example the discussion in |^T[. Because 
of this the goodness of fit depends also on the way of binning the data. To circumvent this 
weakness of the standard goodness of fit test the so-called Parameter Goodness of fit (PG) 



can be used [ffTl. Whereas the standard test measures the probability that all individual 



data points are fitted by an hypothesis, the PG tests the consistency of different data sets 
under an hypotesis. It is based on the function 

XPG = X global, min ,min ; 

(16) 

i 

where Xgiobai.min is minimum of all data sets combined and Xi,min is the minimum of 

the data set i. This function measures the "price" one has to pay by the combination of 
the data sets compared to fitting them independently. It follows a standard distribution 
and should be evaluated for the number of dof corresponding to the number of parameters 
in common to the data sets, see PT| for a precise definition. 



To apply this method we consider the two data sets DAMA versus all the other data 
showing no evidence. Hence, we combine CDMS-Ge, CDMS-Si, and XENON into one data 
set which we denote by NEV. Then we find Xpg ~ 22.6. Evaluating this for 2 dof (corre- 
sponding to the two parameters and ap in common to both data sets) one finds that 
DAMA and NEV data are consistent only at a probability of 1.2 x 10~^. This corresponds 
roughly to the probability of a 2.9a fluctuation in both data sets at the same time.^ We 
conclude that the explanation of DAMA results in terms of spin-independent elastic scatter- 
ing of WIMPs with standard halo properties is strongly disfavoured by XENON and CDMS 
data. Next we investigate the stability of this result with respect to modifications of the 
velocity distribution of the WIMPs in the halo of our galaxy. 



^ Note that the standard Xmin/'iof probabiUty of 5% tests the fit of all individual data points at the best fit 
solution, whereas the PG of 1.2 x 10~^ reflects the compatibility of the DAMA and NEV data sets. These 
are different questions and therefore the two probabilities are consistent with each other. Our DAMA 
analysis is based on 36 data points on the modulation amplitude, where most of them (above 6 keVee) 
are always fitted perfectly, irrespectively of the disagreement with CDMS/XENON. This is an example 
of the above mentioned "dilution" of the standard goodness of fit test. 
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IV. NON-STANDARD HALOS 



The precise limits on the cross section and mass of a dark matter candidate which are 
obtained from a particular observation/non-observation of a signal in a direct detection 
experiment depend upon the velocity dispersion of dark matter around the detector, and 
consequently ultimately upon astrophysical assumptions. The same set of astrophysical 
assumptions are normally made by different experiments so that their results can be com- 
pared with each other, the first of which being that the dark matter halo of the galaxy is 
an isothermal sphere which means a spherically symmetric density distribution of the form 
p oc r~^. For such a density distribution the Keplerian velocity is independent of radius and 
the value of this velocity normally assumed for dark matter studies is a radius independent 
Keplerian velocity of 220 kms~^. It is also assumed that the velocity dispersion of the dark 
matter profile is everywhere isotropic and Gaussian, the width of the Gaussian distribution 
corresponding to the Keplerian velocity of the profile. It turns out that we do not actually 
expect any of these assumptions to hold true for a realistic dark matter halo. 

Over the past decade, A^-body simulations of increasingly large numbers of dark matter 
particles have allowed us to obtain more information about the kind of dark matter halos 



that one would expect to form in an expanding universe (see e.g. |^). These simulations 
show that one might expect a dark matter density that decreases more steeply with radius 
at large radii rather than have the same power law at all radii as in an isothermal profile 
Furthermore, the orbits taken by dark matter particles in a realistic simulation are 
usually rather radial, resulting in an anisotropic velocity dispersion . 



Note that one can also obtain an anisotropic velocity dispersion in a halo where the density 
distribution is not spherically symmetric, but rather triaxial as in ||32|. In this work however, 
we only consider dark matter halos where the density distribution is spherically symmetric 
although the velocity dispersion does not have to be. Also, in a non-extensive ideal gas 
where there is a long range attractive force between the particles such as we have here in 
the form of gravity, one generically expects deviations from Gaussian velocity distribution 
Furthermore, if the dark matter is still not completely virialised but is still coming into 



equilibrium, there will be a superposition of multiple dark matter populations at any given 
place in the halo. This effect would also lead to deviations from a Gaussian distribution of 
velocities. 

In 2006, the results from a Milky Way size dark matter halo simulation called Via Lactea 
containing 234 million particles were published [^|. We have looked at this data to see how 
much the dark matter distribution experienced by an observer at the Solar radius within this 
simulation would vary from the normal assumptions stated above for observers on Earth. 
For each particle in the simulation there is a position vector Xi and velocity vector Vi, plus 
the local gravitational potential per unit mass in units of velocity squared U {xi). We express 
the velocities in terms of the square root of their local potential Vi = Vi/ ^/~U{xi). Next we 
work out the angle between the radial direction and the overall velocity vector. We then use 
this to decompose the velocity into a radial part and a part perpendicular to that which we 
call tangential. For each radius we bin the tangential and the radial velocities obtaining two 
distributions. We fit the one dimensional radial distribution using the following expression 
which we find to be a better fit than the Tsallis distributions (designed to fit non-extensive 
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FIG. 3: The parameters at and fi explained in the text fitted to the radial and tangential velocity 
dispersions of dark matter at different radii from the centre of the Via Lactea simulation. The vertical lines 
indicate the position of the Sun at r = 8.5 kpc. The velocity dispersions are clearly non-Gaussian as one 
approaches the centre of the galaxy. 



or multiple temperature distributions) used in |3 



exp 



R 



3l 



(17) 



Because we have rescaled the velocities with respect to ^J—U{xi), and are dimen- 
sionless constants of order one. The naive assumption is that the width of the velocity 
distribution at a given radius is simply the Keplerian velocity at that radius. While infor- 
mation about the mass distribution of dark matter is required to go from the potential to 
the Keplerian velocity, The parameter //j is an indication of how badly this assumption is 
broken, a encodes the deviation from Gaussianity (a = 1 corresponding to a Gaussian). 
For the one dimensional case, the normalisation is analytic, Nji = 2/Rr(l + l/2a{i). We 
perform the same fitting procedure for the tangential velocity, fitting 



2ttvt 
1^ 



exp 



-2 \ 



(18) 



and while we are not aware of an analytic expression for Nt it is trivial to obtain it nu- 
merically. Note, in terms of the two dimensions perpendicular to the radial direction R, 



o~,2 



r,2 



Using the data from the Via Lactea simulation, we have fitted for values of fi and ctj as a 
function of radius from the centre of the galaxy. The results, which can be seen in figure ^, 
show that there is a considerable deviation from Gaussianity in the velocity dispersion of the 
galaxy. Both the deviation from Gaussianity, the anisotropy of the velocity dispersion and 
the change in the relationship between the width of the dispersion and the local Keplerian 
velocity will change the ratio between the expected modulation in the DAMA experiment 
and the total expected events at XENON and CDMS. We have calculated these changes and 
attempted to see if they can increase the likelihood of the results from both experiments 
being compatible, the results of the new fits can be seen in figure §(a). 



12 




FIG. 4: Allowed regions at 90% and 99.73% CL for DAMA, and exclusion contours for CDMS-Si, CDMS- 
Ge and XENON at 90% CL for the DM halo obtained in the Via Lactea simulation (a), an isotropic 
Maxwellian halo with dispersion w = 110 km/s (b), an asymmetric Maxwellian halo with dispersion vr = 
142 km/s in the radial direction and vt = 63 km/s in the tangential direction (c), and an isotropic Maxwellian 
halo with dispersion v — 220 km/s and escape velocity v^sc — 450 km/s (d). The best fit for DAMA is 
marked with a star. In the panels (b) and (c) we show also the 90% and 99.73% CL regions for the global 
data combining all experiments, as well as the global best fit point (marked with a dot). 

It turns out that the deviation from the assumptions of the isothermal sphere which are 
predicted by the Via Lactea simulation are not sufficient to bring the region in parameter 
space favoured by DAMA away from the region disfavoured by XENON and CDMS. The 
numbers associated with these regions are provided in table H and show that using the 
velocity dispersion predicted by Via Lactea leads to only a very small reduction in and 
the goodness of fit is still unacceptably small. 

It is therefore interesting to ask what kind of halo parameters could lead to a better fit to 
the data, and how realistic would such parameters be? There are examples in the literature 
of the use of a stream of dark matter to boost the annual modulation signal In this 
work, we choose to retain the spherical symmetry of the halo density and instead of adding 
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TABLE I: Summary of the fits to DAMA data and global data (DAMA, CDMS-Ge, CDMS-Si, XENON) 
for different WIMP halos. We give the best fit values, the goodness of fit (assuming 43 dof), the PG 
testing the consistency of DAMA with all other data, as well as the best fit WIMP masses (in GeV) . 

a stream, vary both the width of the velocity dispersion and the anisotropy parameter jS^ei 
defined as 

P,el = l-^. (19) 

A reduction in the width of the velocity dispersion of dark matter alone helps reconcile 
the two data sets without the need to introduce anisotropy. If we assume an isotropic 
distribution of dark matter {Pvei = 0) and reduce v to 110 km/ s which is half of the Keplerian 
velocity at the solar radius, the goodness of fit increases dramatically (see table |ID. A further 
improvement in the fit is made if one assumes a velocity dispersion lower than Keplerian, but 
also highly anisotropic such that Vf( = 142 km/s and vt = 63 km/s. As visible in figure § (c) 
in this cases the entire 90% CL region of DAMA is consistent with the 90% CL bounds from 
CDMS and XENON. For the asymmetric velocity distribution the global drops by about 
20 units compared to the default analysis and provides an excellent goodness of fit of 62%. 
The PG test gives compatibility of DAMA with NEV data with a probability of 4%, due to 
the remaining constraint from XENON. 

A valid question is then whether such low and anisotropic values of the velocity dispersions 
are at all realistic. In order to check on the feasibility of such values, one needs to think about 
particular dark matter halos and see if the solutions of the (Maxwell-) Jeans equations allow 
simultaneously both a high velocity anisotropy (/5^e/ ~ 0.8 ) and low velocity distribution at 
the location of the Sun. 

In order to solve the Jeans equations, we will need to understand the distribution of mass 
in the galaxy. Integration of a spherically symmetric dark matter profile with a well defined 
functional form is trivial. However, at the radius of the Sun, it is important to consider not 
only dark matter but also the presence of baryons, which make up most of the mass in the 
central regions of the galaxy. To model the Milky Way baryon density we assume cylindrical 
symmetry and ignore any spiral arms or bars. For the central bulge of stars we assume a 
density of the form p oc r~"'e~'^^^ while for the disk we assume a (Kuzmin) delta function 
of matter in the z direction {z is the coordinate perpendicular to the disk) with a surface 
density o'diskf'") = '^^^^iskoo^ _ choose the parameters of the model to match observations 

27r(r2+c2)^ 

of the Milky Way: 7 = 1.85, A = Ikpc, c = 5kpc and with the total disk and bulge mass 
Mdiskoo = SMbuige = 6.5 X 1O^°M0 [p], H, 0, H]. We assume that the disk comes to an 
end at a radius of 15 kpc. 

In order to parametrise our dark matter density profile we will consider a profile which 
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assumes two asymptotic radial power law behaviors at both small (7) and large radii^. 
In this profile, known as the 'a/37' profile (or the Zhao profile), the density as a function of 
radius is given by the expression 

P{r) = . . ^ (20) 



(r/a)T [1 + {r/a 



a 



where a governs the radial rate at which the profile interpolates between the asymptotic 
powers —7 and —(3. The parameter a is a characteristic scale radius determining the location 
dividing the two regions described by a single power law. 

Having assumed a value for a, /?, 7 and a we then solve for po in order to get the correct 
value of the Keplerian velocity at the solar radius. This also determines the location of the 
virial radius r^r which in this work is defined to be the radius of the sphere within which 
the average density is 250 times the critical density of the universe (we assume h = 0.7). 
The ratio between r^r and a is referred to as the concentration of the dark matter halo. 

Once we are in possession of these parameters, we can proceed to solve the Jeans equation 
for the radial velocity dispersion 



I djpvl) I 2/3,ez4 _ Vc\ .21) 

p dr r dr r ' 

where 0(r) and Vc{r) are the potential and Keplerian velocity at a given radius. We integrate 
this equation inwards from a large radius several times the magnitude of r^r where we assume 
that pvjj^ = 0. We have checked that the result at r = 8.5 kpc is independent of the exact 
radius at which this boundary condition is applied. We have also assumed that, in the 
absence of a better approximation, the anisotropy parameter /?^e« is a constant throughout 
the halo. 

The results are plotted in figure ^ and show that for a NFW profile where the parameters 
are chosen such that (a,/3, 7) = (1,3, 1) it seems to be rather difficult to imagine that such 
a large value of the velocity anisotropy P^ei could be consistent with low enough values 
of the velocity dispersion to match the data. We also look at a non-standard halo with 
{a, (3, 7) = (1, 4, 1.5). The inner slope of such a halo is quite steep, but even larger values of 
7 than this may be expected in dark matter halos where adiabatic contraction due to the 
presence of baryons has occurred The rate at which density decreases at larger radii is 



also larger than what is normally assumed. 

If we are willing to accept such parameters for the dark matter proffie, it seems that the 
highly anisotropic value of jSyei ~ 0.8 that we require as one ingredient to make the DAM A 
data more consistent with XENON and CDMS is not completely inconsistent with the very 
low velocity dispersions that form the other ingredient. The analysis presented here is meant 
only as a suggestion of the magnitude of possible effects. If such explanations of the DAMA 
data were to be taken seriously, a much deeper analysis of the Jeans equations should be 
undertaken. 

Finally we mention that if one assumes a very low dark matter escape velocity at the solar 
radius then one would remove many of the fastest moving dark matter particles which would 



® This l3 in the density profile should not be confused with the velocity dispersion anisotropy parameter 

Pvel- 
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FIG. 5: Here we plot the radial velocity dispersion vn at the solar radius r = 8.5 kpc as a function of 
the concentration of the dark matter halo Vyir/a for two different dark matter profiles. We have assumed 
that the velocity dispersion anisotropy parameter Py^i = 0.8 and is a constant with respect to radius. The 
horizontal line corresponds to the value of which helps explain the discrepancy. It appears that only for 
sets of halo parameters such as {a,/3,j) — (1,4, 1.5) can one reconcile such a high value of /3 with a low 
enough radial velocity dispersion to help explain the discrepancy between DAMA and XENON/CDMS. 



leave the halo. This would also result in more accord between DAMA and other experiments 
but obtaining a large enough effect is difficult - as expressed in table |, lowering the escape 
velocity to 450 km/s would only marginally make the fit more acceptable. This, however, 
must be considered an unrealistic solution, since the escape velocity at the Solar radius is 
already 440 km/s even if there were no more matter in the Galaxy at larger radii. 

To summarise this section, it seems that without the use of streams but rather by consid- 
ering highly anisotropic velocity dispersions with magnitudes far below the local Keplerian 
velocity at the radius of the Sun would it be possible to reduce the conflict between DAMA 
and XENON/CDMS. 



V. CONCLUSIONS 



Prompted by recent results from DAMA/LIBRA which establish the annual modulation 
of their event rate at the 8.2a level, we have studied the interpretation of this signal in terms 
of spin-independent elastic WIMP scattering. We have shown that the energy spectrum of 
the modulation signal strongly restricts the region of WIMP masses below 10 GeV, confining 
WIMP masses consistent with the DAMA data close to ~ 12 GeV. This region is 
excluded by the limits from CDMS and XENON, and therefore we conclude that even 
if channeling is taken into account this interpretation of the DAMA modulation signal is 
disfavoured. Applying a stringent test to evaluate the consistency of DAMA with null-result 
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experiments we find consistency only with a formal probability of 10~^. 

We have studied how robust this result is with respect to variations of the WIMP velocity 
distribution in our galaxy by changing various parameters of the distribution function. We 
find that decreasing the dispersion of the distribution can somewhat reduce the tension in 
the fit. Adopting in addition an asymmetric WIMP velocity profile with a larger dispersion 
in the radial direction than tangential improves the fit considerably. We conclude that in 
principle it is possible to reconsile DAMA in the considered framework, at the price of rather 
exotic properties of the DM halo. The question remains whether such halo properties can 
be realistic at all. We have checked that a WIMP velocity distribution based on the Via 
Lactea A^-body dark matter simulation does not improve the fit considerably with respect 
to the standard Maxwellian halo model. 

Finally we mention that the negative conclusion on the compatibility of DAMA with 
CDMS and XENON relies crucially on the energy threshold of the latter two. In particular, 
a shift in the nuclear recoil energy scale in these experiments may change the conclusion. 
Indeed, the new measurements of the £efr parameter in XENON |Q (which has not been 
implemented in the first arXiv version of this work) made the disagreement between DAMA 
and XENON somewhat less sever. 
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APPENDIX A: COMMENTS ON THE DAMA SPECTRAL INFORMATION 

Our results are largely based on the fact that DAMA spectral information excludes the 
low-mass WIMP region below 10 GeV. Obviously any effect which affects the spectral shape 
of the signal will have an impact on this conclusion. First, the smearing due to the energy 
resolution of the detector is important. We have checked this by artificially increasing the 
width of the energy resolution function given in Eq. |^ by a factor of two. The global fit 
improves by roughly 7 units in x^, but the tension between DAMA and NEV data persists 
at the level of 8 x 10~^, compare table || and figure ^ (left). 

From figure ^ (left) it follows that the somewhat low data point in the first energy 
bin is very important in constraining the WIMP mass. We have repeated the analysis by 
excluding this bin from the fit, using only the data on the modulation signal above 2.5 keVee. 
In this case the DAMA allowed region extends to lower values of the WIMP mass, and 
once the NEV data are added the globally allowed region includes values of ~ 4 GeV 
and (jp ~ 10~^^ cm^ at 90% CL. This region originates from channeled events on Sodium 
which now can accommodate the spectrum without the first bin, despite the contribution 
of un-channeled Na events. Let us note, however, that in this region constraints from other 



experiments, like CRESST-I fH), TEXONO [|7l, or CoGeNT are relevant. The global 
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TABLE II: Summary of fits to DAMA data and global data (DAMA, CDMS-Ge, CDMS-Si, XENON) 
for the two ad-hoc modifications of the DAMA analysis of figure ^. We give the best fit values, the 
goodness of fit (assuming 42 dof for the last row and 43 dof otherwise), the PG testing the consistency of 
DAMA with all other data, as well as the best fit WIMP masses (in GeV). 




FIG. 6: Allowed regions at 90% and 99.73% CL for DAMA, and exclusion contours for CDMS-Si, CDMS- 
Ge and XENON at 90% CL for two ad-hoc modifications of the DAMA analysis. Left: we artificially assume 
an energy resolution in DAMA a factor two worse than the value given in Right: omitting the lowest 
energy bin of the annual modulation spectrum between 2 and 2.5 keVee. The best fit for DAMA is marked 
with a star. In the right panel we show also the 90% and 99.73% CL regions for the global data combining 
all experiments, as well as the global best fit point (marked with a dot). 



best fit point has an excellent Xgiob min 

= 44.4/42 dof.^ 

Furthermore, we remark that any systematical uncertainty affecting the low energy spec- 
trum may be relevant. For example, figure 26 of shows that the efficiency for DAMA 
single-hit events starts to deviate from 1 below about 8 keVee, just in the signal region. 
Therefore, a possible uncertainty on this low energy efficiency may affect the exclusion of 
the light WIMP region. In the absence of detailed information on possible energy- dependent 
systematic uncertainties we neglect such effects in our analysis. Let us note that the ratio 
of the signals in June and December would be less affected by systematics, since any mul- 



^ The improvement of the consistency of NEV and DAMA data is only partially visible in the PG value 
of about 0.1%, which is still rather low. The reason is that also the fit of DAMA alone improves from 
A^DAMA min ~ '^'^•^ ^'^ iO-i by dropping the first bin, which compensates partially the improvement in the 
global fit. 
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tiplicative uncertainty (even energy dependent) would cancel, whereas the rate difference 
published by the DAMA collaboration is affected by such uncertainties. 

Finally, we mention that the so-called Migdal effect could lead to modifications of the pre- 
dicted energy spectrum in DAMA, see for a discussion and references. An investigation 
of this effect is beyond the scope of this work. 



APPENDIX B: COMPARISON WITH OTHER STUDIES 

In this appendix we compare the results of our work to some studies from other authors. 
The authors of come to a positive conclusion on the consistency between DAMA and 
constraints from other experiments, since they do not include the information on the spectral 
shape of the DAMA signal. In a work which appeared after ours on the preprint server |^ , 



the same authors performed also an analysis including the spectrum which is in agreement 
with our results. 

Our work appeared on the preprint server basically at the same time as [^, with similar 
results. As in our study these authors emphasize the importance of the spectral information 
of the DAMA annual modulation and the constraint from the total unmodulated rate. 
Whereas |^ discusses DM streams, our work considers non-standard halo models. 

A similar study has been performed in stressing also the relevance of spectral in- 
formation and extending the analysis to spin-dependent cross sections. The general results 
for the spin-independent case are in quantitative agreement with us, though in some cases 
the authors draw different conclusions. In particular, they use a variety of statistical tests 
complementary to ours. Whereas our methods are largely based on parameter estimation 
(Ax^ values with respect to the best fit point), these authors show also contours of proba- 
bilities from a goodness-of-fit test based on absolute values. As mentioned at the end of 



section |IIT| this method is often not very sensitive to a tension between different data sets, 
see also footnote ^. Furthermore, by showing contours up to the 5 and even 7a CL they 
do find overlap regions. Let us also comment on the best fit point for DAMA, obtained at 
80 GeV in Tab. IV of [Q, compared to our result = 12 GeV from Eq. |1^. The 



m 



X 



reason why we disfavour the fit in the large DM mass region is the inclusion of the constraint 
from the unmodulated rate in DAMA, which cuts away large part of this region, including 
also the best fit point of [Q, compare figure |l], and shifts the global best fit point to the 
low mass region. 
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